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Location of Impulsive Acoustic Sources in Urban Environments Using Finite- 
Difference, Time-Domain Modeling of Time Reversal with Data from Small Sensor 

Arrays 

Mark E. Todaro, US Army Armament, Research, Development, and Engineering Center, 

Benet Laboratories, Watervliet NY 12189-4000 


ABSTRACT 

The author studied the use of time reversal processing with small sensor arrays, less than ten 
meters across, to locate impulsive sound sources in an urban environment. Sound propagation 
from the source was modeled in two dimensions using a finite-difference, time-domain (FDTD) 
method. The model produced simulated pressure-versus-time traces at sensor locations for 
various sensor arrays. These pressure traces were then reversed in time and fed back into the 
FDTD model to produce a simulated time-reversed wave. In this manner, a number of sensor 
arrays were modeled in various hypothetical urban environments to identify the limitations of 
using time reversal to locate an impulsive source with the smallest arrays possible. Generally, 
arrays of eight or more sensors with a total spatial extent of eight to ten meters were sufficient 
for locating non-line-of-sight sources, whereas smaller arrays proved inadequate. 
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INTRODUCTION AND APPROACH 


Sound propagation in urban terrain is characterized by multiple propagation paths due to 
reflection and diffraction. These effects reduce the reliability of gunfire location systems based 
on beamforming or pulse arrival times, particularly for non-line-of-sight sources. If, however, 
the environment that causes the sound scattering (such as buildings) can be modeled on a 
computer and adequate sensor data can be collected, one can use time-reversal processing to 
model the propagation of a time-reversed waveform returning to the source. Time reversal has 
the added advantage that multiple propagation paths tend to enhance the fidelity of the returning 
wave. 

The preliminary phase of the effort involved creating software to carry out finite-difference, 
time-domain (FDTD) modeling of two-dimensional sound propagation in a region containing a 
number of polygonal reflective structures, such has buildings and walls. The main thrust of the 
effort was then using this software to locate an impulsive sound source in various hypothetical 
urban environments using small sensor arrays and time-reversal processing. The focus was on 
identifying the capabilities and limitations associated with this approach, including the minimum 
workable array size. Particular attention was paid to locating a non-line-of-sight source (around 
the corner or on the other side of a building). 

Time Reversal 


In both acoustics and optics, the wave equation is symmetric with regard to time reversal: 

1 d 2 if/ 


vy- 


dt 


= 0 (where (//is an arbitrary function of space and time, t is time, and c is the 


speed of wave travel). Thus, if y, z,t) is a solution of the wave equation, the time-reversed 
function (//(x, y, z,—t ) is also a solution. The creation of time reversed acoustic waveforms and 
the use of these waveforms to refocus energy at the source is well documented. 1 " 5 


In an idealized demonstration of time reversal, sound would first radiate from a source, as shown 
in Figure 1, propagate through a region containing objects or materials that scatter sound 
(including reflection and diffraction), and be recorded by a number of sensors surrounding the 
region. After the event, the signal from each sensor would be reversed in time and broadcast 
from emitters placed at the corresponding sensor locations. (Conveniently, for many 
applications, piezoelectric ultrasonic transducers can function both as sensors and emitters.) 
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Figure 1. Idealized time reversal. On the left, sound radiates from a source, is scattered, 
and is recorded by closely spaced sensors surrounding the region. On the right, each 
sensor broadcasts a time-reversed signal, creating a time-reversed sound wave that 
returns to the original source. 

According to the Huygens’ principle, the waves from each emitter will superimpose to create a 
waveform identical to the original waveform in all respects except that the propagation vector 
will be at all points directed opposite to the propagation vector of the original wave. Further, 
because of the time symmetry of the wave equation, that reversed waveform will travel (and 
scatter) back through the region and return to the original source. The time-dependent sound 
pressure field throughout the region would be completely reversed as a function of time. Hence, 
the array of sensors combined with the processing to reverse the waveform is sometimes called a 
time-reversal mirror. 

It should be emphasized that either portion of the above process, forward propagation or time- 
reversed propagation, can be made to occur either physically or numerically (computationally 
modeled on a computer). Thus both portions could occur physically, with actual sound traveling 
from a source to the sensors, being time reversed electronically, and then being broadcast from 
the sensors and sent back to the source. In this way, for example, high-energy ultrasound can be 

6 7 

precisely focused to break up a gallstone or kidney stone. ’ 

In the process of locating a sound source in a physical environment, the forward propagation 
would occur physically, the sound would be picked up by sensors, and the time reversed 
propagation would then be modeled on a computer to discover where the reversed wave would 
focus. In this manner, time reversal has been demonstrated for locating gunfire in urban terrain 8 
as well as submarines in the ocean. These type of applications require a pre-existing computer 
model of the acoustic environment. 
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For an ideal time-reversal mirror, the sensors are closely spaced and completely surround the 
area containing the source. For a more practical time reversal mirror, the array of sensors 
subtends a smaller angle as viewed from the source, and we expect the fidelity of the time 
reversed wave to be reduced. If the reduced array were viewed as a traditional mirror, the 
diffraction-limited resolution would then be given approximately by fX! D , where /is the focal 

length (distance to the source in this case), X is the wavelength, and D is the aperture of the 
mirror (transverse size of the array). 

Time reversal, however, makes use of all waves arriving at the sensors, including those arriving 
via multiple propagation paths due to scattering. Thus a wide angular distribution of scatterers 
(as viewed from the source) will direct waves toward the sensors that would otherwise, without 
the widely distributed scatters, not reach the sensors. As a result, multiple scattering typically 
leads to better than diffraction-limited resolution. 10 ’ 11 If the source is effectively immersed in a 
scattering environment, the array of sensors need not surround the environment (as stipulated by 
the ideal), but may be made fairly small. Similarly, one would expect it would be possible to 
reduce the number of sensors within an array of a given size and still maintain fidelity of time 
reversal to an acceptable level. 

The Army’s Cold Regions Research and Engineering Laboratory has studied time reversal to 
locate an impulsive sound source in urban terrain with three to eight sensors distributed 
throughout the environment. They used a 150 meter by 100 meter full-scale artificial urban 
training village with 15 concrete buildings (McKenna MOUT site, Ft. Benning, Georgia). They 
found that even with as few as three sensors, it may be possible to locate the source through time 
reversal processing. It must be noted, however, that although the sensors were few in number, 
they were widely spaced and effectively surrounded, or nearly surrounded, the source 
geometrically. The issue of how small the array (or aperture) might be made was not addressed. 

The effort reported here, however, does address that issue. In this work, time reversal was 
attempted for relatively small sensor arrays, ten meters across or less, completely remote from 
the source rather than “enclosing” it in any way. As viewed from the source, then, the array 
aperture was relatively small. 

Finite-Difference, Time-Domain Modeling 


It is more straightforward to apply the finite-difference method to acoustics by starting with the 
continuity and Euler’s equations rather than the wave equation. 

In two dimensions, the linear continuity equation can be written as 


du dv 1 dp 
- 1 - 1 --- 

ax dy pc dt 



where u is the particle velocity in the x-direction, v is the particle velocity in the y-direction, p is 

the acoustic pressure (instantaneous pressure minus the constant equilibrium pressure in the 

12 

medium), p is the constant equilibrium density, and c is the speed of sound in the medium. 
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The two components of the linear Euler’s equation (also called the force equation) are given by 


12 


dp du dp dv 

— + p — = 0 and — + p — = 0. 


dx dt 


dy dt 


These equations can then be written as the following difference equations: 


A p = -pc 2 At 


r Au Av 2 
- 1 - 

V Ax Ay ) 


(discrete continuity equation) 


A u = - 


At 

P 


r Ap^ 

v Ax J 


(x component of discrete Euler’s equation) 


Av = - 


At 

P 


( Ap^ 

\ A yj 


(y component of discrete Euler’s equation) 


For speed of computation, the finite-difference calculation is carried out using the so-called 
staggered Yee mesh. In this approach, pressures are calculated at discrete points on a square 
lattice and the particle velocities at points on shifted lattices, u shifted half a lattice vector in x 
and v shifted half a lattice vector in y, as shown in Figure 2. 



-©-©-©-©- 



-(p)—(u)—(p)—(u)—(p y- 

Figure 2. Two-dimensional staggered Yee mesh used in the finite-difference 
calculation. Pressures are calculated at lattice points labeled p. The particle velocities 
are calculated at points labeled u and v. 


In addition, the pressure and particle velocity are calculated in alternating time steps. Thus, in 
one time step, the pressures are calculated from the adjacent velocity components (available 
from the previous step) via the discrete continuity equation. At the next time step, the velocity 
components u and v are calculated from the adjacent pressures via the discrete Euler’s equation. 
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Reflections due to truncation of the domain at its outer edges are avoided by use of an absorbing, 
nonreflecting boundary based on Yuan et al.’s method 14 of implementing Berenger’s perfectly 
matched layer. 15 The perfectly matched layer in this case has eight layers. 

In finite difference processing, reflection, refraction, and diffraction occur as a result of spatial 
variations in the equilibrium density p and the speed of sound c. In this application, reflection 
and diffraction of sound at solid barriers such as walls and buildings is of primary concern, and 
the transmission of sound within such objects are ignored. Hence a vast simplification will be 
realized by assuming these solids do not support shear waves and are thus “fluid” in the acoustic 
sense. The change of acoustic impedance {pc) at the air-solid interfaces will be sufficiently high 
that a negligible amount of transmitted acoustic energy re-emerges from the solids. 

Because the high contrast of acoustic impedance at the air-solid interface can result in solution 
instability, one-lattice density averaging is used at those interfaces following Schroder and 
Scott. 16 


Reflection and diffraction at an interface (between two fluids, as discussed above) is determined 
by the change in acoustic impedance pc at the interface. For simplicity, the acoustic impedance 
for the structures (walls and buildings) in these simulations will be assumed to be that of 
concrete, a highly reflective material with both a higher density and a higher sound velocity than 
air. For finite difference modeling, however, high sound velocities are undesirable due to the so- 
called Courant condition for stability. For a square lattice, the time step At must satisfy 


A t< 



to guarantee numerical stability (where c max is the highest sound velocity in the model). 17 (This 
condition can be intuitively understood as stating that during time-stepping, information must 
travel faster than sound.) A higher sound velocity in any region of the model thus requires a 
shorter time step (or larger lattice spacing). It is not the sound velocity per se that determines 
reflection and diffraction, however, but the product pc. Hence the desired acoustic impedance 
for the structures can be achieved by keeping the sound velocity equal to that of air and 
compensating by assigning an artifically high density. 


Thus 344 m/s, the sound velocity of dry air at 21° C, is used for both air and the structures. A 
density of 1.2 kg/m 3 is used for air. The density of the structures is assumed to be 19,440 kg/m 3 , 
to give a structure-to-air impedance ratio of 16,200, similar to what would be given by the actual 
velocity and density of concrete (2950 m/s and 2300 kg/m 3 ). 
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For the impulsive source, the model uses a point source of pressure given by the following 
waveform, which has been shown to reasonably approximate the pressure pulse from a gunshot 
or explosion: 8 ’ 18 


Pit) = A 




exp (-Cft) 


where A is an arbitrary amplitude, B = 0.25, C=3.39, and/= 150 Hz. To speed the calculation 
slightly, P is assumed to be zero after 0.025 s. This waveform is plotted in Figure 3. 



Figure 3. Waveform used for impulsive sound source. 


For the finite-difference calculations, the spacing between adjacent pressure points (and adjacent 
velocity points) was 0.1 m, less than 1/20 of the nominal wavelength of the source function. The 
calculations used a time step of 0.1 ms between successive pressure calculations (and between 
successive velocity calculations). 

The calculations were carried out using MATLAB. For evaluation and visualization purposes, 
each run of the finite-difference calculations produced a video file (AVI) of the pressure field 
evolution and a file containing pressure versus time for various sensor positions. 
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RESULTS 


Figure 4 shows the calculated acoustic pulse propagating from a point source in one hypothetical 
urban area that was studied. Although other hypothetical regions were studied, this region is 
convenient for illustrating a number of typical results. In all cases shown for this region, each 
calculation was carried out for 350 ms of sound travel to ensure the complete waveform and 
reflections were received at the sensors. 



Figure 4. Calculated acoustic pulse propagating from a point source in a 92 meter by 68 
meter hypothetical urban area. The four images are snapshots of the sound pressure at 
different times. Buildings are white and open areas are gray. Positive sound pressure is 
depicted as white and negative pressure as black. The source is shown as a single white 
circle at the bottom of the vertical “street.” A seven-sensor hexagonal array eight meters 
across is shown as seven white circles in the horizontal “street.” 
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The images in Figure 4 also show a seven-sensor hexagonal array eight meters across (two 
meters between adjacent sensors). The pressure traces were recorded at these sensor positions 
for the subsequent time reversal processing shown in Figure 5. The two lower images of Figure 
5 show the time-reversed wave converging at the location of the original source, focusing at the 
bright white spot in the lower right image. Note that time reversal with this hexagonal array was 
able to locate a non-line-of-sight source. This array, the largest studied, generally produced 
good results. 



Figure 5. Time-reversed sound from seven-sensor hexagonal array. The two lower 
snapshots show the time-reversed wave converging at the location of the original source, 
focusing at the bright white spot in the lower right image. (Source as in Figure 4.) 


8 









Figure 6 shows four snapshots of the time-reversed pressure field for the same point source as in 
Figure 4 but with an array of two small subarrays ten meters apart. Each subarray consisted of 
four sensors with one at the center and the other three each one meter from the central sensor and 
arranged on an equilateral triangle. Again, the lower right image shows the time-reversed sound 
focusing at the original source location. 



Figure 6. Time-reversed sound from an array of two small subarrays ten meters apart. 
The two lower snapshots show the time-reversed wave converging at the location of the 
original source, focusing at the bright white spot in the lower right image. (Source as in 
Figure 4.) 
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Figure 7 shows four snapshots of the time-reversed pressure field for the same point source as in 
Figure 4 but with a five-sensor portion of the hexagonal array close to the wall of a building 
which acts as a mirror for sound. 



Figure 7. Time-reversed sound from a five-sensor portion of the seven-sensor hexagonal 
array backed against a wall. The two lower snapshots show the time-reversed wave 
converging at the location of the original source, focusing at the bright white spot in the 
lower right image. (Source as in Figure 4.) 
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Figure 8 shows the calculated acoustic pulse propagating from a point source at a different 
location from previous examples. 



Figure 8. Calculated acoustic pulse propagating from a point source at a different 
location from previous examples. 


Figures 9 and 10 show the time-reversed pressure field for two different sensor arrays using the 
sound from the source at the second location. 



Figure 9. Time-reversed sound from the seven-sensor hexagonal array. The two 
snapshots show the time-reversed wave converging at the location of the source. 
(Source as in Figure 8.) 
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Figure 10. Time reversed sound from the array of two small subarrays ten meters apart. 
(Source as in Figure 8.) 


Figure 11 shows the calculated acoustic pulse propagating from a point source at a third location. 



Figure 11 .Calculated acoustic pulse propagating from a point source at a third location. 


Figures 12 and 13 show the time-reversed pressure field for two different sensor arrays using the 
sound from the source at the third location. 



Figure 12. Time-reversed sound from the seven-sensor hexagonal array. 
(Source as in Figure 11.) 
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Figure 13. Time-reversed sound from the array of two small subarrays ten meters apart. 
(Source as in Figure 11.) 


The previous examples have all been of sensor arrays locating non-line-of-sight sources. In 
locating non-line-of-sight sources, however, there must be a reflective geometry that provides 
reasonably strong sound propagation from the source to the sensors. Hence an array on one side 
of a building would have difficulty locating a source on the opposite side unless there were 
additional buildings beyond the first to reflect sound toward the sensor area. 


The next two figures illustrate the location of a direct, line-of-sight source using the seven-sensor 
hexagonal array. Figure 14 shows the calculated acoustic pulse propagating from the point 
source located at the left of the region. 



Figure 14. Calculated acoustic pulse propagating from a line-of-sight point 
source. 
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Figure 15 then shows snapshots of the time-reversed wave approaching and converging at the 
original source. 



Figure 15. Time-reversed sound from the seven-sensor hexagonal array converging to a 
line-of-sight source. (Source as in Figure 14.) 


Figures 16 and 17 illustrate how multiple propagation paths enhance time reversal. The left side 
of Figure 16 shows a sequence of snapshots of the calculated acoustic pulse propagating from a 
point source with no objects to scatter sound. The right side shows the calculated acoustic pulse 
propagating from an identical point source but with a reflecting corridor formed by two parallel 
walls. In both cases, a linear array of six sensors is located at the right of the region, 40 meters 
from the source and 2 meters between adjacent sensors. The walls are 29 meters long and 16 
meters apart. 
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elapsed time = 30.0 ms elapsed time = 30.0 ms 





elapsed time = 120.0 ms 



elapsed time = 120.0 ms 



Figure 16. Snapshots of the calculated acoustic pulse propagating from a point source 
with no scattering objects (on the left) and with a reflective corridor (on the right). A 
linear array of six sensors is at the right of both regions. 


Figure 17 shows the time-reversed waves from the six-sensor array with and without the walls. 
In the presence of the walls, the sound refocuses much more tightly at the original source 
location. 
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Figure 17. Time-reversed waves from the six-sensor array with and without 
the walls. The sound refocuses much more tightly with the walls present. 
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As stated earlier, the array of sensors combined with time-reversal processing is sometimes 
called a time-reversal mirror. If the reduced array is viewed as a traditional mirror, the 
diffraction-limited resolution is approximately fA/ D , where /is the distance to the source, A is 

the wavelength, and D is the transverse size of the array. The corridor walls, however, direct 
sound toward the sensors that would otherwise be lost. Hence, the walls effectively increase the 
aperture D of the time-reversal mirror, thus reducing the size of the refocused spot. 

Consistent with that interpretation, it was generally observed that adding buildings and reflective 
propagation paths tended to enhance the ability of time reversal to locate the source. 

CONCLUSION 

The use of time reversal processing with small sensor arrays to locate impulsive sound sources 
was studied by calculating the forward and time-reversed propagation of sound in various 
hypothetical urban environments. The focus was on relatively small sensor arrays, ten meters 
across or smaller, completely remote from the source rather than “enclosing” it in any way, a 
scenario far from ideal for time reversal. Although this leads to poor fidelity of the overall time- 
reversed wave, it is still often possible to determine the location of the source from a video 
animation. 

Generally speaking, the results confirmed that the ability to locate a source decreases as the 
number of sensors and the spatial extent of the array is reduced. Part of the purpose of the 
investigation was to learn just how small an array could be made yet still be able to locate a 
source, particularly non-line-of-sight. Arrays of seven or more sensors with a total spatial extent 
of eight to ten meters were sufficient for locating non-line-of-sight sources, whereas smaller 
arrays proved inadequate. A seven-sensor hexagonal array eight meters across, for example, 
performs quite well in locating a source in urban terrain, even on the opposite side of a building. 
A more practical array consisting of two smaller four-sensor subarrays ten meters apart is still 
adequate for non-line-of-sight and around-the-corner source location. (Because the distance 
between subarrays is not critical, so long as it is approximately ten meters or more, each subarray 
can conceivably be carried on a vehicle.) 

As the size of the arrays is reduced below eight to ten meters, the performance degrades 
significantly. A single small array of three or four sensors spaced one to two meters apart, for 
example, is generally insufficient to locate a source. 

One practical drawback of the FDTD method is its long processing time, even in two 
dimensions. For a typical 90x70 square meter area modeled in this study, a 3.6 GHz personal 
computer took 15 minutes to process 350 milliseconds of sound propagation. 

Given the high degree of reflection in an urban environment, however, it may be possible to 
neglect diffraction when using time reversal to locate the source. Hence a less computationally 
intense approach based on geometrical acoustics, such as beam tracing, 19 ’ 20 may be feasible. A 
beam tracing approach could make time reversal practical for source location (1) with small 
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sensor arrays as discussed in this report; (2) with synchronized sensors distributed throughout a 
town (already demonstrated by others using time reversal with FDTD processing); and (3) with 
three-dimensional modeling for rooftop or tower-based sensor arrays. 

Although the results presented here are not from actual physical measurements, they do illustrate 
the potential that time reversal holds for source location. 
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